High-quality metagenome-assembled genomes from proximal colonic microbiomes of synbiotic-treated korean native black pigs reveal changes in functional capacity

Synbiotics are feed supplements with the potential to promote health and productivity in pigs partly, through modulation of the intestinal microbiome. Our study used shotgun sequencing and 16S rRNA gene sequencing techniques to characterize the effect of a synbiotic containing three Lactobacillus species and a fructo-oligosaccharide on the proximal colonic microbiome of 4- to 7-month-old Korean native black gilts. With shotgun sequencing we constructed unique metagenome-assembled genomes of gut microbiota in Native Black Pig for the first time, which we then used for downstream analysis. Results showed that synbiotic treatment did not alter microbial diversity and evenness within the proximal colons, but altered composition of some members of the Lactobacillaceae, Enterococcaceae and Streptococcaceae families. Functional analysis of the shotgun sequence data revealed 8 clusters of orthologous groups (COGs) that were differentially represented in the proximal colonic microbiomes of synbiotic-treated Jeju black pigs relative to controls. In conclusion, our results show that administering this synbiotic causes changes in the functional capacity of the proximal colonic microbiome of the Korean native black pig. This study improves our understanding of the potential impact of synbiotics on the colonic microbiome of Korean native black pigs.

The Jeju Black Pig is a small, black domestic pig breed that is native to the Korean island of Jeju where it has been conserved for generations 1 . The breed is well adapted to foraging in the rough landscape of the island, where it endures windy conditions and an overall harsh climate 2 . Meanwhile the evolutionary pressure in modern commercial pig farming systems on the mainland has led to the dominance of exotic and commercially developed breeds. Modern commercial breeds have been bred for high performance traits such as prolificacy, mothering ability, number of piglets born, weight gain, growth rate, back fat thickness, lean meat, among others. Despite a low genetic potential for weight gain and overall performance in Jeju black pigs, their pork products have a unique taste which makes them highly sought after and relatively pricier on the market 3,4 . In recent years, there have been efforts to transform Jeju pork production from a traditional backyard system to a more productive, industrialized system while continuing to preserve the breed. Innovative management strategies are being employed to optimize the productivity of this native pig breed. Among these strategies is the use of probiotics, prebiotics and synbiotics to promote health and weight gain in the farmed pigs.
Synbiotics are a combination of probiotics and prebiotics that are used as feed additives or supplements in order to grant health benefits to the host 5 . Probiotics have been defined as "live microorganisms that, when administered in adequate amounts, confer a health benefit on the host" 6 . Prebiotics are dietary components, usually plant fiber that are undigestible by the host but selectively support the growth or activity of a limited group of beneficial microbial genera or species in the gut 7 . These synbiotic combinations modify intestinal microbial composition and influence the development and functioning of the gut with putative benefits to health and performance of pigs. Probiotics containing members of the genus Lactobacillus have been used to modulate the intestinal microbial flora and promote gut health, weight gain and overall performance of pigs [8][9][10] . This genus is a core member of the pig's intestinal microbial ecosystem 11 . Lactobacilli are of interest because of their inhibitory role on some potentially pathogenic bacteria in the swine gut. They produce lactic acid which has a disruptive effect on the walls of gram-negative bacteria, such as the E. coli and salmonella which are key pathogens in domestic pigs 12,13 . This in addition to their production of antimicrobial peptides 14 and modulation of immune response in pigs 15 makes lactobacillus-based probiotics potentially good alternatives to antibiotics in pig farming. Besides contributing to the host's immune response to enteric pathogens, some strains of Lactobacillus promote nutrient digestibility by producing dietary enzymes such as phytases and amylases which then enables breakdown of indigestible dietary components 16 . Several studies have shown that combining prebiotics such as fructo-oligosaccharides with probiotics synergistically enhances the benefits of the probiotic 17,18 . These effects are most important in the colonic environment.
As a site with the highest number of microorganisms, the proximal colon plays very important roles not only in the excretion of waste but also in the absorption of water and electrolytes as well as absorption of vitamins produced by its complex and rich microbial flora 19 . A healthy balance of the microbial ecosystem within the colon is central to its role in the swine gut in that colonic mucosa is the main surface where the gut microbiota and the immune system interact 20 . Therefore, manipulation of this environment to promote the health of Native Black pigs is of interest. The potential of probiotics to promote the productivity of these pigs, can be understood further, by studying their effects on the proximal colonic microbial ecosystem. Metataxonomic techniques such as targeted sequencing of the 16S ribosomal RNA gene and metagenomic techniques like, shotgun sequencing, are frequently used to study microbial composition and infer biological functions in such communities 21,22 . The former approach is cheaper since, by targeting one gene (and typically, a fragment of it), the method imposes lower demands on sequencing coverage, compared to the latter, which probes all prevailing genes in an environment. However, techniques that target marker genes, are limited when studying functional profiles of microbial communities 22 . Even though computational techniques can be used to predict functionality from say, 16SrRNA gene-based taxonomies 23 , shotgun sequencing offers better resolution of functional gene composition within microbial environments 22 . Use of shotgun sequencing data is highly dependent on the completeness of the reference databases for the environment in question. To date, however, the microbiome of Jeju black pig has not been sufficiently studied to support a metagenomic assembly approach.
In this study, we used both shotgun and 16S rRNA sequencing techniques to explore, for the first time, the proximal colonic microbiota of Jeju Native Black pigs. We used a bioinformatic approach to construct a custom database for the Jeju Black pig proximal colonic samples to improve taxonomic classification of shotgun sequence data. We then evaluated the effect of a synbiotic formulation containing Lactobacillus buchneri NLRI-1201, L. plantarum NLRI-101, L. casei DK128 and the prebiotic, fructo-oligosaccharide on the microbial composition within the proximal colon.

Methods
Animals and sample collection. All animal experiments were conducted in a facility of contract research organization at Cronex Co., Ltd. (Cheongju, Korea), with following the regulations of the Cronex Co., Ltd. Animal Care and Use Committee (CRONEX-IACUC 202002007). Jeju black pigs (M-Pig®) were purchased from Cronex Co., Ltd. and housed in Cronex Animal SPF (Specific Pathogen Free) facility for whole experimental period, where was environmentally controlled with temperature and relative humidity at 22 ± 2 °C and 50 ± 10%, respectively, under a 12-h light-dark cycle. Each cage (4 pigs per cage, 1.0 m 2 per pig) was equipped with plastic mesh floor, two round feeders, and two nipple drinkers. Samples were collected from 7-month-old female Jeju 16S rRNA metagenome sequencing. From the isolated gDNA, the V3-V4 hypervariable region of the 16S rRNA gene was amplified using the universal prokaryotic primers Bakt_341F (5′-CCT ACG GGNGGC WGC AG-3′) and Bakt_805R (5′-GAC TAC HVGGG TAT CTA ATC C-3′) 24 . Amplicon libraries were prepared, pooled, and loaded onto an Illumina MiSeq flow-cell using a 600-cycle reagent cartridge (v3 reagent kit). Sequencing was done using a 2 × 301 bp paired-end sequencing method. 16S rRNA metagenome analysis. Denoising and merging, the read pairs were processed using DADA2 package (q2-dada2 v2020.8.0) 25 (with options -p-trunc-len-f 275 -p-trunc-len-r 225). Denoised features were classified by Qiime2 q2-feature-classifier plugin 26 using a naïve Bayes classifier trained on the V3-V4 hypervariable region of 16S rRNA gene sequences in the Silva Database v138 27 . This denoising and taxonomic classification process was done using the Qiime2 (v2020.08) software 28 . Alpha and beta diversity metrics were calculated using normalized/rarefied sub-samples (14,397 features) of the denoised reads. Analysis of alpha and beta diversity was achieved using Qiime2 diversity core-metrics-phylogenetic plugin and results were visualized in R 29 using the package ggplot2 30 . Following Thorsen and colleagues' 31 recommendations for small datasets, we normalized the features using Trimmed Means of M-values (TMM), and explored differentially abundant taxonomy using edgeR 32 .
Host genome sequences were removed by extracting pair-unmapped reads. To reconstruct metagenome-assembled genomes (MAGs) from our samples, the extracted reads were used as input for de novo assembly using the metaSPAdes v3.14.1 tool 35 .

Bioinformatic analysis of read-based taxonomies.
A read-based taxonomic classification was conducted using kraken2 v2.1.1 software 41 . We used both the in-built Standard Kraken2 database, and a custom host database to classify the paired end reads from the Jeju Black pig colonic microbial environment. The custom host database was built based on the default database and included the National Center for Biotechnology Information's (NCBI) RefSeq genomes as well as our classified MAGs. For comparison, reads were aligned to Choc-oPhlAn database (2019.01v) using MetaPhlAn3 (v3.0.4), and UniRef90 database using HuManN3 (v3.0.0) 43 .
Following classification in Kraken, and with a normalized read length parameter set at 150, the relative taxonomic abundances of microbes in the samples were re-estimated using bracken v2.6.2 software 44 . Alpha and beta diversity analyses were calculated and visualized using the R package, phyloSeq (v1.34.0) 45 . Differential abundance of the microbes in the samples was statistically analyzed using the edgeR package 32 in R. Only species with a Log-transformed count-per-million value > 1 and a False Discovery rate (FDR) < 0.01 were considered differentially abundant with statistical significance.
Gene-based metagenome assembly. For gene-based differential abundance analysis, PLASS assembly (release 4-687d7) 46 was employed to obtain protein sequences for each sample. Assembled protein sequences were merged, dereplicated and clustered using the easy-linlinclust option in MMseqs2 (version 13.45111) 47 with sequence identity parameter set at 0.5. Sequence searches were then performed using DIAMOND 2.0.11 software 48 and with the option "blastx -evalue 0.001 -fast -max-target-seqs 1", gene cluster abundances were estimated. Clustered representative sequences were then functionally annotated using emapper-2.1.5 49 basing on eggnog orthology data 50 . The R package, edgeR 32 was then used to estimate group specific changes in relative

Results
16S rRNA metagenome sequencing. To investigate possible difference in taxonomic composition between the groups, V3-V4 hypervariable regions of the 16S rRNA gene were sequenced from microbial communities in the proximal colon of the study animals. Paired-end sequences were quality-filtered, merged and denoised leaving > 10,000 high quality reads. We assessed possible difference in alpha diversity between the two groups based on 16S rRNA metagenome sequences. Alpha diversity analysis using Shannon and Simpson indices indicated that there was no significance in within-group diversity between our symbiotic-fed animal (treatment group) and controls. The mean number of Shannon index in the treatment group was 10.447 ± 0.135, and 9.896 ± 0.903 for control group. Also, Simpson index was 0.9989 ± 0.0002 and 0.9982 ± 0.0011 for treatment group and control group, respectively (Fig. 1A). Analysis of beta diversity distances using Bray-Curtis dissimilarity showed that symbiotic-fed group is well clustered while control group showed sparse distribution (Fig. 1B). Taxonomic analysis indicated that Firmicutes and Bacteroidetes were the most abundant phyla found in the samples. On average, firmicutes had a relatively lower abundance in the treatment group (62.8% ± 9.288) compared to the control group (74.1% ± 3.764). Meanwhile, members of the phylum, Bacteroidetes had a higher relative abundance (25.3% ± 3.815) in the treatment group compared to controls (18.6% ± 4.728) (Fig. 1C). However, www.nature.com/scientificreports/ differential abundance in genus level testing between groups found no statistical significance (FDR p-value < 0.05) in the observed differences.
Shotgun metagenome sequencing. Sequencing on MGISEQ generated 116 Giga Bases with mean quality scores of 34.8 from our samples. After quality trimming and removal of host-derived sequences, about 12.7 Giga Bases (a mean of 43 million read pairs) per sample remained. Single-sample assembly and binning of resulting contigs produced 1,603 bins, which were then quality checked in CheckM (Fig. S2). Concatenation, filtration, and dereplication of all the bins resulted in a total of 360 high-quality metagenome-assembled genomes (MAGs). Taxonomies assigned by a GTDB-based classification of the MAGs as well as the phylogenetic relationships between the MAGs are presented in Fig. 2A Read-based classification of the shotgun sequences based on the Standard Kraken2 Database yielded a low average alignment (25.07%) (Fig. 2B). To improve classification rate, we built a custom host database which put together our 360 MAGs and the Standard Kraken DB. With this, we were able to achieve a higher average read alignment rate (84.49%). Also, Metaphlan3/Humann3 read classification analysis were conducted and compared with Kraken2 results, showing that almost all of reads were not aligned in nucleotide alignment and about 50 ~ 60% of reads were not aligned in translated alignment (Table S1).
With the read-based taxonomies, we estimated the alpha-and beta-diversities of the microbiota within the proximal colons of the Jeju black pig gilts. By estimating colonic microbial diversities using our shotgun dataset, we obtained lower alpha diversities in the treatment group relative to controls (Fig. S3A). This is in contrast to earlier results from diversity analysis using the 16S rRNA dataset which showed the reverse; higher alpha diversities in the treatment group compared to the controls group (Fig. 1A). However, data from both sequencing methods showed larger deviations in the alpha diversity values within the treatment group compared to the control group. Analysis of Beta diversity distances found no group specific clustering (Fig. S3B).

Differential abundance analysis.
Considering the indistinguishability in abundance and diversity between groups at community level, further analysis was conducted to explore differential abundance at the spe-  www.nature.com/scientificreports/ cies level (Fig. 2C). In total, 6 species were identified as statistically different between the treatment group and the control group. Among these species, 5 were bacteria belonging to the firmicutes phylum, while the other was a bacteriophage (Streptococcus virus sfi21) (Fig. 2D) Functional profiling of pig microbiome. To study the antimicrobial resistome of the proximal colonic samples from our Jeju black pig gilts, we scanned for antibiotic resistance genes within the metagenome assembled contigs. From the samples, we found a total of 73 antibiotic resistance genes (Fig. 3A). The distribution of these AMR genes varied among individuals and showed no differential abundances between groups. We then investigated the functional role and diversity of proteins encoded by microbial genes associated with colonic microbiomes of the pig samples. To get insights into the gene cluster abundances in metagenomes, gene-based assembly was conducted, and 90,555,907 non-redundant coding sequences were generated. We then assigned functional categories to the genes by mapping them to Clusters of Orthologous Groups (COG) categories. Overall, results showed very little variation in the number of COGs between groups, and samples had coterminous COG profiles (Table S2). Results also showed that 7 COG categories of known functions had average abundances above 5.0% within the samples and together accounted for 53.56% of the total gene abundances. These categories included the COG L (replication & repair), COG G (carbohydrate metabolism and transport), COG M (Cell wall/membrane/envelop biogenesis), COG E (amino acid metabolism and transport), COG J (Translation), COG K (Transcription) and COG C (Energy production and conversion). Differential analysis uncovered 12 differentially abundant genes from 8 COG groups that passed the statistical significance cut-off of FDR corrected p-values < 0.0005 (Table 2, Fig. 3B).
Gene-based assignment of taxonomies yielded phylum-level profiles presented in (Fig. 3C). The proportion of the phylum, Bacteroidetes was 20.37% (± 0.016), and that of Firmicutes was 53.18% (± 0.017). These proportions were in agreement with the 16S rRNA dataset. Protein-level profiles were nearly the same between samples.
We also mapped our data to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database to yield KEGG orthologues. Analysis of the KEGG orthologs revealed a statistically significant increase in the 'Iron complex transport system substrate-binding protein' within the treatment group. While the relative increase of the 'mucolipin 3' and '3-aminobutyryl-CoA ammonia-lyase' groups in the control group was statistically significant. We however, found no statistically significant differences based on other functional categories (CAZyme, KEGG pathway and PFAM).

Discussion
In this study, we looked to improve our understanding of the effect of a synbiotic on the colonic microbiome of Korean native black pigs that are traditionally reared on the Jeju Island. We combined shotgun metagenomic sequencing and 16s rRNA marker gene sequencing techniques to study the impact on the microbial composition in the proximal colon of these pigs. The results indicate that the synbiotic causes a subtle shift in the functional profile of the proximal colonic microbiome but no detectable shift in phylogenetic diversity. This is also the first study into the colonic microbiome of this native breed of pigs that are raised under the unique environmental conditions of Jeju Island.
During analysis of the shotgun metagenome sequence data, we encountered a low rate of alignment in the existing reference databases and tools (Fig. 2B, Table S1). One of the reason for this could be the limitation in coverage of a database that is still under curation, considering the quantity of data required to meticulous record all complete bacterial, archaeal, and viral genomes in all environments 52 . Secondly, we suspect that due to the spatial www.nature.com/scientificreports/ isolation of this native breed of pigs on Jeju Island, their gut environments might represent unique, previously unstudied microbiomes which are not included in the kraken database. To overcome this limitation, we built a custom database for our data through identifying novel MAGs. Consequently 360 novel and high-quality MAGs were identified from the proximal colonic luminal content samples in this study ( Fig. 2A). Classification based on the customized database resulted in high-resolution taxonomic composition, with fewer unclassified reads. Although 10 ~ 20% of the reads remained unclassified, we were able to explain more of the unknown reads from the sequence data. In this way, the analysis managed overcome the limitations imposed by a reference database that could not reflect many of the rarely researched microbiota of Native Black Pig. We found that dietary supplementation with a synbiotic formula had no statistically significant effect on the microbial evenness, alpha-and beta-diversity in the proximal colon of Korean native pigs. The age of the piglets might be one of the reasons why the synbiotic did not have a detectable impact on the community-level microbial diversity 8 . The developing gut of young piglets is more susceptible to significant shifts and becomes less so as it develops into a richer, more mature and stable environment. Our study animals were weaner-growers or sub-adults which have a relatively stable gut microbial composition compared to younger piglets in the periweaning phases of growth. Although there was no obvious clustering according to study group, we found some unique/ differentially abundant taxa between the groups. There was an increased representation of some genera of the Lactobacillaceae family (Weissela spp., Pediococcus spp.,) and Enterococcaceae (Enterococcus spp.) probably due to the stimulatory effect of the probiotic Lactobacilli 53 or the fructo-oligosaccharides in the synbiotic supplement 54 . There was also a reduction in some species of the Streptococcaceae family (Streptococcus spp.) as well as a streptococcus targeting bacteriophage (Streptococcus virus sfi21).
In the functional profiling of the proximal colonic microbiomes, we used both shotgun sequence and 16S rRNA marker gene sequence data (Fig. 3). We were able to detect some statistically significant effects of the synbiotic treatment on the distribution of microbial functional gene groups using shotgun sequence data but not with the 16S rRNA sequence data. This was not surprising, as it has been reported in several previous studies that the shotgun sequence data yields more genes and functional metagenomes than the 16S rRNA sequence data 55,56 . During functional profiling, the shotgun sequence data enabled identification of antimicrobial resistance (AMR) genes in the proximal colonic microbiomes of our study animals. We found a total of 73 AMR genes in the proximal colonic microbiome of the study animals, and their abundance did not exhibit group-specific shifts in this study (Fig. 3A). Genes responsible for resistance against tetracyclines, macrolides, aminoglycosides, lincosamides and beta lactams were found in the majority if not all of the samples. This finding in Jeju Black pigs studied here, is similar to findings from several studies that have reported occurrence of resistance genes among commercial pig breeds 57,58 . Pollock and colleagues 57 , found that occurrence and abundance of these AMR genes, however, did not always correspond to antimicrobial exposure. In our study, the animals were not exposed to antibiotics throughout the experimental period, however exposure prior to inclusion into the study was not investigated. Nonetheless, AMR genes are known to be ubiquitous among microbial communities, having coevolved with bacterial strains as a strategy to counter toxic chemicals within their niches 59 .
In addition to identifying antimicrobial resistance genes within this microbial community, the metagenomic analysis provided insights into synbiotic-associated shifts in the functional capacity of the proximal colonic microbiomes (Fig. 3B). The analysis revealed 12 statistically significant differentially distributed gene (belonging to 8 COG groups) between the study groups. Synbiotic treated pigs were enriched in gene categories involved in secondary metabolites biosynthesis, transport, and catabolism (COG Q), Signal transduction mechanisms (COG T) and Replication, recombination, and repair (COG L). The COG Q category includes genes encoding Isochorismatase, an ether hydrolase that catalyzes the synthesis of iron-binding molecules called siderophores. Siderophores have been reported to act as virulence factors 60,61 , and in some instances to mediate microbial interactions that result in suppression of potential gut pathogens 62,63 . Genes in the COG T enable microorganism to sense, and respond to a number of intracellular and extracellular stimuli 64 . Therefore, an abundance of gene families involved in signal transduction could indicate a higher capacity for the bacteria in this group to respond to changes in their environment including nutrient changes, antibiotics, pH among others 65 . While the COG L family of genes is involved in DNA transposition in bacteria. An abundance of this gene family might indicate a high potential for genetic recombination 66 and/or lateral gene transfer in the microbial community within the proximal colon of synbiotic-fed pigs 67 . Interestingly, some of the relatively increased genes in the treatment group, originated from Lachnospiraceae bacterium and Dorea longicatena, which are stimulated by fructo-oligosaccharides 68,69 .
We also found that gene families involved in coenzyme transport and metabolism (COG H), transcription (COG K), translation, ribosomal structure, and biogenesis (COG J) as well as carbohydrate transport and metabolism (COG G) were significantly lower in the synbiotic-fed pigs compared to controls. Genes in the COG H category are involved in coenzyme transport and metabolism. Coenzymes are involved in a wide range of intracellular reactions. They play important roles in energy metabolism, cellular response to stress, intracellular redox homeostasis, nucleic acid repair, xenobiotic metabolism, and sometimes, modulation of microbial virulence [70][71][72][73][74][75] . Genes belonging to the COG K category code for RNA polymerase sigma factors and are involved in initiation of bacterial transcription during exponential growth and gene transcription in response to environmental stressors as well as during morphological changes 76,77 . Fewer COG K genes in the proximal colonic microbiome of synbiotic treated pigs might imply a relative reduction in bacteria with the capacity for rapid multiplication and appropriate response to new environmental stimuli. The COG J category of genes codes for the ribosomal L22 protein which is a structural component of the 50S ribosomal subunit where it binds with multiple domains in the 23S rRNA and facilitates ribosomal assembly 78 . The L22 protein forms part of the peptide exit tunnel of bacterial ribosomes and is therefore thought to regulate translation 79 . The protein also interacts with and mediates susceptibility to macrolides which work by inhibiting translation 80 . COG G is a family of genes that encodes extracellular solute-binding proteins which facilitate bacterial uptake and utilization of carbohydrates 81 www.nature.com/scientificreports/ lower representation could mean a reduced ability to utilize carbohydrates in their environment including the prebiotic studied here 83 . It is worth noting that the some of the significantly increased genes within the control groups, originated from alkaliphilic bacteria (Indibacter alkaliphilus, Halanaerobium hydrogenifromans, Paenibacillus daejeonensis). A possible reason for this might be that relative increase of these genes in the control group microbiomes could be a result of an actual reduction of these genes in the microbiomes of the treatment group and not necessarily an increase in the controls. This is due to the acidic conditions created by the lactic acid producing bacteria in the synbiotic and the relative proliferation of lactobacillus producing bacteria that is stimulated in the treatment group. Overall, the mechanistic association between the synbiotic and the significant differentially expressed genes between the groups is unclear. However, the pattern of differential distribution suggests a shift of functional capacity away from stress tolerance within the synbiotic fed group. This reduced capacity/need for stress tolerance, in turn, could indicate a stabilizing effect of the synbiont on the environment conditions within the proximal colons of the treated pigs.
A major limitation of this study was the small sample size within study groups. Nonetheless, we hope that when interpreted within these limits, the study provides insights into this previously unexplored microbial environment of an unstudied Korean native breed of pigs.

Conclusion
In this study, we used metagenomic techniques to study the impact of a synbiotic on the proximal colonic microbiomes of Korean native, Jeju black pig gilts. Our results show that the synbiotic, containing Lactobacillus buchneri NLRI-1201, L. plantarum NLRI-101, L. casei DK128 & fructo-oligosaccharide does not cause significant shifts in the microbial abundances, evenness, and diversity at community level. However, the synbiotic increases abundance of members of the family, Lactobacillaceae which have been found to offer health promoting effects to their hosts. Further functional analysis uncovered an effect of the synbiotic on genes belonging to 12 COGs within the colonic microbiome of the Jeju black pigs implying changes in functional capacity of the colonic microbiome. This work improves our understanding of the effect of synbiotics on the microbiome of rare pig breeds such as the Korean black pigs that are native to Jeju Island. We can also infer an understanding of the synbiotics' potential to promote productivity and industrialization of Korean black pigs.